function [ A_monot ] = A_monot_fun(combos_all_unique,J)

%% sigma is increasing in: own x, own z
%%          decreasing in: other x other z

n_Bcoeff = size(combos_all_unique,1);

mat = cell(n_Bcoeff,1);
sub_minus = cell(n_Bcoeff,1);

for c=1:n_Bcoeff

    temp_vec = combos_all_unique(c,:);
    temp_geq = all(repmat(temp_vec([1 J+1]),[size(combos_all_unique,1),1])<=combos_all_unique(:,[1 J+1]),2);             % take combos with HIGHER own x, own z        
    temp_leq = all(repmat(temp_vec([2:J J+2:2*J]),[size(combos_all_unique,1),1])>=combos_all_unique(:,[2:J J+2:2*J]),2);     % take combos with LOWER other x, other z      
    temp_final = logical(temp_geq.*temp_leq);
    sub_minus{c} = find(temp_final);
    sub_minus{c} = setdiff(sub_minus{c},c);

end

for c=1:n_Bcoeff
    
    sub_minus_union = [];
    
    for i=1:length(sub_minus{c})
        
        sub_minus_union = union(sub_minus_union,sub_minus{sub_minus{c}(i)});
        
    end

    temp_sub_col_minus = setdiff(sub_minus{c},sub_minus_union);             % eliminate redundant inequalities, i.e. lines (because of transitivity)
    temp_sub_row_minus = [1:length(temp_sub_col_minus)]';
    temp_mat = zeros(length(temp_sub_col_minus),n_Bcoeff);
    
    if isempty(temp_sub_col_minus)==0
        
    temp_ind_minus = sub2ind(size(temp_mat),temp_sub_row_minus,temp_sub_col_minus);
    
    temp_mat(temp_ind_minus) = -1;
    
    temp_sub_col_plus = c*ones(length(temp_sub_col_minus),1);
    temp_sub_row_plus = [1:length(temp_sub_col_minus)]';
    
    temp_ind_plus = sub2ind(size(temp_mat),temp_sub_row_plus,temp_sub_col_plus);
    
    temp_mat(temp_ind_plus) = +1;
    
    temp_mat((temp_sub_col_minus==c),:) = [];     % eliminate row for c
    
    end
    
    mat{c} = sparse(temp_mat);
    %n_Bcoeff-c
end

A_monot = cell2mat(mat);

A_monot = unique(A_monot,'rows');

end